Week 6 Computer Lab Worksheet

Multilevel models

Author

Chris Moreh

Session aims

  • learn using the lmer() function from the lme4 package to fit multilevel models with two or more levels of nesting
  • dealing with survey weights in R

Exercise 1: Economic and Cultural Drivers of Immigrant Support

In this exercise we will replicate a model presented in Table 4 of Valentino et al. (2019) (first column, titled “ALL”), and extend its multilevel framework. Valentino et al. (2019) explore the causal determinants of public opposition to immigration in 11 countries using a series of survey treatments relative to type of employment, racial composition, family status, and country of origin. The data consists of repeated measures taken from the same survey respondents (subjects), and the authors employ “mixed models” (their term) to account for multiple observations per subject, as these are likely to be correlated given that they belong to the same respondent.

So, we have clustering at two levels: at the individual level and at a country level. The models fitted by the authors take into considerations and adjust for the individual-level clustering, but not the clustering at the country level. In the model that we will first replicate, the authors pool together respondents from all the countries (see the models named ALL in the regression tables), while also including separate regression models for each country. The “ALL” model allows the identification of global trends that persist regardless of country-specific features; however, it does not account for the fact that respondents - and therefore the errors - are also affected by country-level influences. The second part of the exercise consists of fitting a multilevel model that accounts for both clustering factors. The multilevel model contains information that summarises both the “ALL” and the various separately run country-level regression models.

1.1. Two-level random intercepts model

First, we will load the original dataset and perform a number of data-cleansing tasks to set up the data in a similar format to that used by the authors:

# Load required packages
pacman::p_load(tidyverse, lme4, sjlabelled)

# Import the dataset
val <- sjlabelled::read_stata("https://cgmoreh.github.io/HSS8005/data/Valentino17.dta")
# Filter out cases
val <- val |> 
  filter(!is.na(uid)) |>      # remove rows with no uid
  filter(!is.na(support))     # remove rows where dependent variable is not provided

# Create country codes to match to country name key
countryt_key <- data.frame(countryt = c("AU", "CA", "CH", "DK", "ES", "FR", "JP",
                                        "KR", "NO", "UK", "US"),
                           country_name = c("Australia", "Canada", "Switzerland",
                                            "Denmark", "Spain", "France", "Japan", 
                                            "S.Korea", "Norway", "UK", "USA"))

# Creating new variables
val <- val |> 
  mutate(candf = factor(cand),     # create factor versions of the treatment variables  
         Job_f = factor(treat_hstat, levels = c(0, 1), labels = get_labels(treat_hstat)),
         Complexion_f = factor(treat_drk, levels = c(0, 1), labels = get_labels(treat_drk)),
         ME_f = factor(treat_nat_me, levels = c(0, 1), labels = get_labels(treat_nat_me)),
         Children_f = factor(treat_kids, levels = c(0, 1), labels = get_labels(treat_kids)),
         Gender_f = factor(treat_gender, levels = c(0, 1), labels = get_labels(treat_gender)))

# Recode gender variable to ease interpretation, and add country name
val <- val |> 
  mutate(female = factor(gender, levels = c(1,0), labels = c("female", "male"))) |> 
  left_join(countryt_key, by = "countryt")

# Set other categorical variables as factors

## It's useful to first make a list of the variable names:
fact_vars <- c("id", 'uid', 'country', "countryt", 'cand', 'insample', 
               'gender', 'education2', 'occupation2', "inc2", 'inc3',
               'treat_hstat', 'treat_drk', 'treat_nat_me', 'treat_kids', 'treat_gender')

#  Then recode those listed variables as factor variables; 
# and we can also filter the data to treat_gender equals 0 to follow what the authors did; this should leave us with 37,579 cases/respondents in the dataset, the total N that appears in Table 4, Model "ALL"
val_mod <- val |> 
  mutate_at(fact_vars, factor) |> 
  filter(treat_gender == 0)

To inspect and better understand the dataset, apply some of the summary functions we practised in previous weeks to this data (e.g. using modelsummary, gtsummary, summarytools or something else).

...

We are now ready to fit a multilevel model with random/varying intercepts for each respondent to account for intra-class correlation within individuals:

# Varying-intercept model with an individual-level predictors
## Same model as ALL in Table 4, but with REML = T
m_all <- lmer(support ~ Job_f + Children_f + ME_f + Complexion_f + candf + 
                (1|uid),                         # This is where we add the clustering variable coding individual ID numbers
                 data = val_mod, 
                 weights = wt_country, REML = T) # Some additional settings to emulate results from models fit in Stata

1.2. Fully nested model: repeated obs nested within countries

We can now add an additional level of clustering, including countryt as an additional level:

m_ct <- lmer(support ~ Job_f + Children_f + ME_f + Complexion_f + candf + 
                    (1 | countryt / uid),         # This is where we add additional levels
                  data = val_mod, 
                  weights = wt_country, REML = T) # Some additional settings to emulate results from models fit in Stata

We cna now tabulate the results and compare them. modelsummary::modelsummary() is again a good option, but sjPlot::tab_model() includes some additional transformations to the random effects components out of the box that can be rather useful for a quick glance (we can have more control over those computations and calculate them and include them manually in another table format):

# Modelsummary table
modelsummary::modelsummary(list(m_all, m_ct), stars = TRUE)
 (1)   (2)
(Intercept) 0.523*** 0.528***
(0.004) (0.021)
Job_fHigh status 0.120*** 0.120***
(0.004) (0.004)
Children_fChildren −0.002 −0.002
(0.002) (0.002)
ME_fME −0.020*** −0.020***
(0.002) (0.002)
Complexion_fDark treatment 0.001 0.002
(0.004) (0.004)
candf2 −0.038*** −0.038***
(0.002) (0.002)
SD (Observations) 0.137 0.137
SD (Intercept uid) 0.287
SD (Intercept uidcountryt) 0.279
SD (Intercept countryt) 0.067
Num.Obs. 37579 37579
R2 Marg. 0.039 0.039
R2 Cond. 0.822 0.822
AIC 7392.9 6443.2
BIC 7461.1 6520.0
ICC 0.8 0.8
RMSE 0.11 0.11
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
 
# SJplot table
sjPlot::tab_model(m_all, m_ct, digits = 3)
  DV:Support DV:Support
Predictors Estimates CI p Estimates CI p
(Intercept) 0.523 0.515 – 0.532 <0.001 0.528 0.488 – 0.569 <0.001
Job f: High status 0.120 0.112 – 0.129 <0.001 0.120 0.112 – 0.129 <0.001
Children f: Children -0.002 -0.006 – 0.002 0.253 -0.002 -0.006 – 0.002 0.255
ME f: ME -0.020 -0.023 – -0.017 <0.001 -0.020 -0.023 – -0.017 <0.001
Complexion f: Dark
treatment		 0.001 -0.007 – 0.010 0.748 0.002 -0.006 – 0.011 0.592
candf: candf 2 -0.038 -0.041 – -0.035 <0.001 -0.038 -0.041 – -0.035 <0.001
Random Effects
σ2 0.02 0.02
τ00 0.08 uid 0.08 uid:countryt
  0.00 countryt
ICC 0.81 0.81
N 19734 uid 19734 uid
  11 countryt
Observations 37579 37579
Marginal R2 / Conditional R2 0.039 / 0.822 0.039 / 0.822

Exercise 2: Return to Osterman

At the end of Worksheet 3 we fitted the full Model 1 reported in Table 3 of Österman (2021) (see also Table A.3. in their Appendix for a more complete reporting on the models). Comparing the model we fitted to the one reported by the author, our results were very close, but not totally equivalent on several coefficients and standard errors. The main reasons for the divergence had to do with two aspects of the published model that we had disregarded: (1) we did not include survey weights to correct for sampling errors, and (2) we did not allow for intra-group correlation in the standard errors among respondents from the same country (and the same age cohort). We will first implement these two additional steps and compare our final results again to those reported in Österman (2021). Then, we will refit the model in a multilevel framework.

Refitting Model 1

As a quick reminder, let’s refit the model as we have done on Worksheet 3.

To be noted

There was a mistake in the code used to fit the model on Worksheet 3. This has now been corrected there and below.

In the footnotes to Table 3, the author tells us that in addition to what is reported in the table,

All models include reform FEs [Fixed Effects], a general quadratic birth year trend, and reform-specific trends for birth year (linear) and age (quadratic). Additional controls: foreign-born parents and ESS round dummies.

The version of the model summary presented in Table A.3 in the Appendix also includes the additional controls, but not the other terms.

In the main text, the author further explains some of the reasoning behind their modelling choices:

One dilemma for the design is that there has been a trend of increasing educational attainment throughout the studied time period, which means that the reform-windows of treated and non-treated cohorts will also pick up the effects of this trend. To counter this, [the list of covariates] includes a general quadratic birth year trend, reform-specific linear controls for birth year and reform-specific quadratic age trends. The quadratic terms are included to allow sufficient flexibility in absorbing possible non-linear trends of increasing education within the reform-window of seven treated and seven untreated birth year cohorts. … The reform-fixed effects are also essential because they imply that only the within-reform-window variation is used to estimate the effects and between-reform differences are factored out, such as pre-reform differences in social trust. (Österman 2021, 221–22)

Before we fit the model, some of the concepts in the quotation need unpacking. A quadratic term is a second-degree polynomial term - put simply, it’s the square of the variable concerned. The quadratic of a variable such as \(age\) is therefore \(age^2\), or \(age \times age\). In other words, it is like a variable’s “interaction” with itself. Because of this, there are several ways in which the quadratic terms to be included in a model can be specified:

  1. We could create the quadratic terms as new variables, and include those in the model, as below:
# First, load the original dataset
osterman <- sjlabelled::read_stata("https://cgmoreh.github.io/HSS8005/data/osterman.dta")

# Create new quadratic variables by multiplication with themselves and add them to the dataset, saving it as a new data object
osterman2 <- osterman |> 
  mutate(agea_quad = agea*agea,        # quadratic age variable
         yrbrn_quad = yrbrn*yrbrn)     # quadratic birth-year variable

# We now have two additional variables in the dataset; we can fit the model using those:
m1_prequad <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d +                         # main covariates reported in Table 3
                   fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) +     # additional controls for foreign-born parents and ESS Round dummies
                   agea + yrbrn + agea_quad + yrbrn_quad +                              # general quadratic birth year trend and quadratic age
                   factor(reform_id_num) +                                              # reform fixed effects dummies
                   yrbrn:factor(reform_id_num) +                                        # reform-specific birth year trend
                   agea:factor(reform_id_num) +  agea_quad:factor(reform_id_num),       # reform-specific quadratic age trend
               data = osterman2)                                                        # the new expanded dataset
  1. We can get the same results by creating the quadratic terms directly as part of the modelling function. The one thing we should keep in mind is that if we want to include standard mathematical operations within a formula function, we need to isolate or insulate the operation from R’s formula parsing code using the I() function, which returns the contents of I(...) “as.is”. The model formula would then be:
m1_funquad <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d +                         # main covariates reported in Table 3
                 fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) +       # additional controls for foreign-born parents and ESS Round dummies
                 agea + yrbrn + I(agea^2) + I(yrbrn^2) +                                # general quadratic birth year trend and quadratic age
                 factor(reform_id_num) +                                                # reform fixed effects dummies
                 yrbrn:factor(reform_id_num) +                                          # reform-specific birth year trend
                 agea:factor(reform_id_num) +  I(agea^2):factor(reform_id_num),         # reform-specific quadratic age trend
              data = osterman)                                                          # the original dataset
  1. In the two previous options, the quadratic terms will be correlated with the original variables. To avoid this by relying on so-called orthogonal polynomials we should use the poly() function. We can also fit the same correlated polynomial model as the ones above by adding the raw = TRUE option to the poly() function. In the code below, we will fit the correlated version first, then the orthogonal version (This stackoverflow discussion explains in more detail the difference between the two options):
m1_polyraw <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d + 
                 fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + 
                 poly(agea, 2, raw = TRUE) + poly(yrbrn, 2, raw = TRUE) +
                 factor(reform_id_num) +           
                 yrbrn:factor(reform_id_num) + 
                 agea:factor(reform_id_num) + poly(agea, 2, raw = TRUE):factor(reform_id_num),
               data = osterman)

m1_polyorth <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d + 
                 fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + 
                 poly(agea, 2) + poly(yrbrn, 2) +
                 factor(reform_id_num) +           
                 yrbrn:factor(reform_id_num) + 
                 agea:factor(reform_id_num) + poly(agea, 2):factor(reform_id_num),
               data = osterman)

It’s worth noting, however, that the Stata routine used by the author fitted correlated/raw coded polynomials, so the orthogonal version below is just for a comparison and we will not use it going forward. For a side-by-side overview comparison of the models we fitted so far, we could use a model summary tabulation function (e.g. sjPlot::tab_model() or modelsummary::modelsummary(), which we are already familiar with from previous lab exercises; some other popular options include stargazer::stargazer() and jtools::export_summs(); you can read through the function documentation and test them out). Below we will use modelsummary():

# It's cleaner to first make a list of the models we want to summarise; we can even name them:
models <- list(
  "Pre-calculated quadratic" = m1_prequad,
  "Within-function quadratic" = m1_funquad,
  "poly() with raw coding" = m1_polyraw,
  "poly() with default orthogonal coding" = m1_polyorth
)

# modelsummary table with stars for p-values added
modelsummary::modelsummary(models, stars = TRUE)

The results from the modelsummary() are not shown here because it’s a long and ugly table, but it’s useful for perusing to compare the results across the different models. For a cleaner table showing only the results included in Table A.3 in the Appendix to Österman (2021), we can use the coef_map or the coef_omit option in modelsummary() and only include m1_funquad, which will be the polynomial fitting routine that we will use going forward:

# It's cleaner to first make a vector of the coefficients we wish to include; we can name the coefficients as they appear in Table A.3; note that we also leave out the Intercept, as in the published table:
coefs <- c("reform1_7" = "Reform",
           "female" = "Female",
           "blgetmg_d" = "Ethnic minority",
           "fbrneur" = "Foreign-born father, Europe",
           "mbrneur" = "Foreign-born mother, Europe",
           "fnotbrneur" = "Foreign-born father, outside Europe",
           "mnotbrneur" = "Foreign-born mother, outside Europe",
           "factor(essround)2" = "ESS Round 2",
           "factor(essround)3" = "ESS Round 3",
           "factor(essround)4" = "ESS Round 4",
           "factor(essround)5" = "ESS Round 5",
           "factor(essround)6" = "ESS Round 6",
           "factor(essround)7" = "ESS Round 7",
           "factor(essround)8" = "ESS Round 8",
           "factor(essround)9" = "ESS Round 9")

# Then we pass the vector to coef_map to select the coefficients to print
modelsummary::modelsummary(list(m1_funquad), stars = TRUE, coef_map = coefs)
 (1)
Reform 0.063*
(0.027)
Female 0.058***
(0.013)
Ethnic minority −0.241***
(0.054)
Foreign-born father, Europe −0.111**
(0.042)
Foreign-born mother, Europe −0.108*
(0.044)
Foreign-born father, outside Europe −0.065
(0.073)
Foreign-born mother, outside Europe −0.110
(0.078)
ESS Round 2 0.059
(0.045)
ESS Round 3 0.162*
(0.075)
ESS Round 4 0.243*
(0.108)
ESS Round 5 0.360*
(0.144)
ESS Round 6 0.397*
(0.179)
ESS Round 7 0.449*
(0.212)
ESS Round 8 0.655**
(0.246)
ESS Round 9 0.816**
(0.283)
Num.Obs. 68796
R2 0.200
R2 Adj. 0.198
AIC 268913.8
BIC 270056.1
Log.Lik. −134331.877
RMSE 1.71
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

We can compare our table to the published one to confirm that the results are similar, but there are some inconsistencies that we can address.

Survey weights

One of the differences between our coefficients and the ones in the published table is due to the fact that we did not account for survey weights, while in Österman (2021), p. 220:

The data are weighted using ESS design weights

To understand in more detail what this means, we need some understanding of how weights are constructed in the European Social Survey (ESS). In a nutshell, ESS data are distributed containing 3 types of weights:

  1. dweight: These are the so-called design weights. Quoting the ESS website: “the main purpose of the design weights is to correct for the fact that in some countries respondents have different probabilities to be part of the sample due to the sampling design used.” These weights have been built to correct for the coverage error, that is, the error created by the different chances that individuals from the target population are covered in the sample frame.

  2. pspwght: These are the post-stratification weights. According the the ESS website, these “are a more sophisticated weighting strategy that uses auxiliary information to reduce the sampling error and potential non-response bias.” These weights have been computed after the data has been collected, to correct from differences between population frequencies observed in the sample and the “true” population frequencies (i.e. those provided by the Labour Force Survey funded by the EU and available on Eurostat). Unlike the design weights, which are based on the probability of inclusion of different groups of individuals in the sample frames, these have been calculated starting from variables that are there in the data, and are really an “adjustment” of the design weight made to reach observed frequencies that match those of the target population.

  3. pweight: These are the population size weights. These weights have the purpose to match the numbers of observations collected in each country to the country populations. They are to be used only when we calculate statistics on multiple countries (for instance, unemployment in Scandinavia). Their value is the same for all observations within the same country.

There is a lot to be said about survey weights and options for dealing with them, which we will not cover in more detail in this course. But as part of this exercise we will get to know some functions that can help with including survey weights and can be extended to include more complex design weights as well. Specifically, we will look at the survey package:

## Install and load required packages

options(repos = list(CRAN = "http://cran.rstudio.com/")) # Just in case

pacman::p_load(tidyverse, survey)     

# For some reason installing/loading the tidyverse installed an earlier version of the dplyr package and "srvyr" complained; this installs the latest version of "dplyr":
# install.packages("dplyr")

We start by creating a weighted data object using the svydesign() function from survey, which includes the dweight as used in the original article:

## Create weighted data
osterman_w <- svydesign(id = ~1,                 # specifying cluster IDs is needed; ~0 or ~1 means no clusters
                        weights = ~dweight,      # we apply the design weights
                        data = osterman)

We then fit the model using the svyglm() function from the survey package and save the model object; note that we specify a design = option with the weighted data object we created earlier:

m1_w <- svyglm(trustindex3 ~ reform1_7 + female + blgetmg_d + 
           fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + 
           agea + yrbrn + I(agea^2) + I(yrbrn^2) +
           factor(reform_id_num) +           
           yrbrn:factor(reform_id_num) + agea:factor(reform_id_num) +  
           I(agea^2):factor(reform_id_num),
       design = osterman_w, data = osterman)

To compare the results from the weighted model to the one we produced earlier, we can check them in a modelsummary() table; we can reuse the list of coefficients to display that we created earlier:

# List and name the models
models <- list(
  "Unweighted model" = m1_funquad,
  "Weighted model" = m1_w)

# Tabulate the models
modelsummary::modelsummary(models, stars = TRUE, coef_map = coefs)
Unweighted model Weighted model
Reform 0.063* 0.063*
(0.027) (0.029)
Female 0.058*** 0.061***
(0.013) (0.014)
Ethnic minority −0.241*** −0.261***
(0.054) (0.066)
Foreign-born father, Europe −0.111** −0.090*
(0.042) (0.046)
Foreign-born mother, Europe −0.108* −0.092*
(0.044) (0.046)
Foreign-born father, outside Europe −0.065 −0.053
(0.073) (0.079)
Foreign-born mother, outside Europe −0.110 −0.087
(0.078) (0.082)
ESS Round 2 0.059 0.052
(0.045) (0.047)
ESS Round 3 0.162* 0.148+
(0.075) (0.080)
ESS Round 4 0.243* 0.253*
(0.108) (0.114)
ESS Round 5 0.360* 0.364*
(0.144) (0.153)
ESS Round 6 0.397* 0.403*
(0.179) (0.190)
ESS Round 7 0.449* 0.458*
(0.212) (0.224)
ESS Round 8 0.655** 0.671*
(0.246) (0.261)
ESS Round 9 0.816** 0.835**
(0.283) (0.299)
Num.Obs. 68796 68796
R2 0.200 0.206
R2 Adj. 0.198 0.204
AIC 268913.8
BIC 270056.1 283455.7
Log.Lik. −134331.877 −141031.683
RMSE 1.71 1.71
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

We can check our results again to those published in Table A.3 in the Appendix to Österman (2021), and we will see that all the coefficients are the same now. The remaining differences are in the standard error estimates, and those are due to the fact that we did not use robust standard errors to account for intra-group correlation.

Variance estimation using robust clustered errors

Österman (2021), p. 222 write that:

All models are estimated with OLS, … and apply country-by-birth year clustered robust standard errors. [Footnote: An alternative would be to cluster the standard errors on the country level but such an approach would risk to lead to biased standard errors because of too few clusters]

Applying clustered robust standard errors is a more elementary and less flexible way to account for breaches of the iid assumptions of OLS (that our variables are “independent and identically distributed”) than fitting a mixed-effects (multilevel, hierarchical) model. There is a vast literature assessing whether one approach is better suited than the other in different contexts. To get a feel for the issues in question and for a deeper understanding of how these methods are useful, compare the analysis of Cheah (2009), whose results “suggest that modeling the clustering of the data using a multilevel methods is a better approach than fixing the standard errors of the OLS estimate”, to that of McNeish, Stapleton, and Silverman (2017), who discuss a number of cases (focusing on psychology literature) where cluster-robust standard error may be more advantageous than multilevel/hierarchical modelling.

We will first fit our model using clustered robust standard errors, as done by Österman (2021), then we will check the results against those we would obtain from a multilevel model design to see how the risk related to the low number of country clusters affects the results.

The standard was of applying error corrections is by using the sandwich and {lmtest} packages, but they do not handle well the summary objects produced by the survey package for weighted estimates. Using the unweighed model we fitted earlier, we could do the following:

# Load required packages
pacman::p_load(sandwich, lmtest)

# extract variance-covariance matrix with clustered correction
vc_cl <- vcovCL(m1_funquad, type='HC1', cluster=~cntry_cohort)

# get coefficients
m1_funquad_cl <- coeftest(m1_funquad, vc_cl)

## Or, the two steps above can be combined into one call:

m1_funquad_cl2 <- coeftest(m1_funquad, 
                          vcovCL, type='HC1', cluster = ~cntry_cohort)

# And we can check that the two have the same result:
all.equal(m1_funquad_cl, m1_funquad_cl2)
[1] TRUE

Let’s tabulate the results for comparison:

# List and name the models
models <- list(
  "Unweighted model" = m1_funquad,
  "Weighted model" = m1_w,
  "Unweighted model with cluster-robust errors" = m1_funquad_cl)

# Tabulate the models
modelsummary::modelsummary(models, stars = TRUE, coef_map = coefs)
Unweighted model Weighted model Unweighted model with cluster-robust errors
Reform 0.063* 0.063* 0.063*
(0.027) (0.029) (0.029)
Female 0.058*** 0.061*** 0.058***
(0.013) (0.014) (0.015)
Ethnic minority −0.241*** −0.261*** −0.241***
(0.054) (0.066) (0.060)
Foreign-born father, Europe −0.111** −0.090* −0.111*
(0.042) (0.046) (0.046)
Foreign-born mother, Europe −0.108* −0.092* −0.108*
(0.044) (0.046) (0.047)
Foreign-born father, outside Europe −0.065 −0.053 −0.065
(0.073) (0.079) (0.077)
Foreign-born mother, outside Europe −0.110 −0.087 −0.110
(0.078) (0.082) (0.086)
ESS Round 2 0.059 0.052 0.059
(0.045) (0.047) (0.044)
ESS Round 3 0.162* 0.148+ 0.162*
(0.075) (0.080) (0.079)
ESS Round 4 0.243* 0.253* 0.243*
(0.108) (0.114) (0.115)
ESS Round 5 0.360* 0.364* 0.360*
(0.144) (0.153) (0.154)
ESS Round 6 0.397* 0.403* 0.397*
(0.179) (0.190) (0.190)
ESS Round 7 0.449* 0.458* 0.449*
(0.212) (0.224) (0.224)
ESS Round 8 0.655** 0.671* 0.655*
(0.246) (0.261) (0.262)
ESS Round 9 0.816** 0.835** 0.816**
(0.283) (0.299) (0.298)
Num.Obs. 68796 68796 68796
R2 0.200 0.206
R2 Adj. 0.198 0.204
AIC 268913.8 406007.8
BIC 270056.1 283455.7 1033594.4
Log.Lik. −134331.877 −141031.683
RMSE 1.71 1.71
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

The error-correction does not have an impact on the coefficients, but it did affect the standard errors. Notice that the standard errors came closer to those in the weighted model without any error correction, showing that the weighting procedure already applies an error correction, albeit not specifically on the cluster variable cntry_cohort.

An easy way to include information on the clustering directly in the svydesign() function is to add cluster IDs:

# Re-specify the survey design
osterman_w_cl <- svydesign(id = ~cntry_cohort,        # This time we add cntry_cohort as id
                         weights = ~dweight, 
                         data = osterman)

# Re-fit the model with the new survey design
m1_w_cl <- svyglm(trustindex3 ~ reform1_7 + female + blgetmg_d + 
               fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + 
               agea + yrbrn + I(agea^2) + I(yrbrn^2) +
               factor(reform_id_num) +           
               yrbrn:factor(reform_id_num) + agea:factor(reform_id_num) +  
               I(agea^2):factor(reform_id_num),
           design = osterman_w_cl, data = osterman) 

# Tabulate the models

# Add the latest model to the existing list using the append() function to save typing
models <- append(models, 
                 list("Weighted with cluster IDs" = m1_w_cl)
                 )

# Tabulate the models
modelsummary::modelsummary(models, stars = TRUE, coef_map = coefs)
Unweighted model Weighted model Unweighted model with cluster-robust errors  Weighted with cluster IDs
Reform 0.063* 0.063* 0.063* 0.063*
(0.027) (0.029) (0.029) (0.030)
Female 0.058*** 0.061*** 0.058*** 0.061***
(0.013) (0.014) (0.015) (0.015)
Ethnic minority −0.241*** −0.261*** −0.241*** −0.261***
(0.054) (0.066) (0.060) (0.067)
Foreign-born father, Europe −0.111** −0.090* −0.111* −0.090+
(0.042) (0.046) (0.046) (0.047)
Foreign-born mother, Europe −0.108* −0.092* −0.108* −0.092+
(0.044) (0.046) (0.047) (0.047)
Foreign-born father, outside Europe −0.065 −0.053 −0.065 −0.053
(0.073) (0.079) (0.077) (0.078)
Foreign-born mother, outside Europe −0.110 −0.087 −0.110 −0.087
(0.078) (0.082) (0.086) (0.081)
ESS Round 2 0.059 0.052 0.059 0.052
(0.045) (0.047) (0.044) (0.045)
ESS Round 3 0.162* 0.148+ 0.162* 0.148+
(0.075) (0.080) (0.079) (0.088)
ESS Round 4 0.243* 0.253* 0.243* 0.253*
(0.108) (0.114) (0.115) (0.125)
ESS Round 5 0.360* 0.364* 0.360* 0.364*
(0.144) (0.153) (0.154) (0.169)
ESS Round 6 0.397* 0.403* 0.397* 0.403+
(0.179) (0.190) (0.190) (0.209)
ESS Round 7 0.449* 0.458* 0.449* 0.458+
(0.212) (0.224) (0.224) (0.246)
ESS Round 8 0.655** 0.671* 0.655* 0.671*
(0.246) (0.261) (0.262) (0.288)
ESS Round 9 0.816** 0.835** 0.816** 0.835*
(0.283) (0.299) (0.298) (0.326)
Num.Obs. 68796 68796 68796 68796
R2 0.200 0.206 0.206
R2 Adj. 0.198 0.204 −209.204
AIC 268913.8 406007.8
BIC 270056.1 283455.7 1033594.4 20460875.2
Log.Lik. −134331.877 −141031.683 −10229741.419
RMSE 1.71 1.71 1.71
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

This final model takes us the closest to the model reported in Österman (2021).

To take that a final step further, we will fit the model in a multilevel framework using the lme4 package. In the next exercises, you will rely on your knowledge of the lmer() function to fit the required models, then add them to the modelsummary() table and compare them to the OLS models we fitted earlier.

Exercise 3: Fit a random intercept model

Complete the code below in order to fit a random intercept model with cntry_cohort as the single grouping variable and disregarding survey weights:

# Load the required packages
pacman::p_load(lme4)

# Fit a random intercept model with `cntry_cohort` as the grouping variable
m2 <- ...


# Fit a random intercept model with only `cntry` as the grouping variable
m2_cntry <- ...


# Add the latest models to the existing list using the append() function
models <- append(...)

# Tabulate and compare the models
modelsummary::modelsummary(...)

Questions

  • How do the models compare?

References

Cheah, B. 2009. “Clustering Standard Errors or Modeling Multilevel Data ?” In.
David, F. N. 1955. “Studies in the History of Probability and Statistics i. Dicing and Gaming (a Note on the History of Probability).” Biometrika 42 (1/2): 1–15. https://doi.org/10.2307/2333419.
El-Shagi, Makram, and Alexander Jung. 2015. “Have Minutes Helped Markets to Predict the MPC’s Monetary Policy Decisions?” European Journal of Political Economy 39 (September): 222–34. https://doi.org/10.1016/j.ejpoleco.2015.05.004.
Gelman, Andrew, Jennifer Hill, and Aki Vehtari. 2020. Regression and other stories. Cambridge: Cambridge University Press. https://doi.org/10.1017/9781139161879.
Lord, R. D. 1958. “Studies in the History of Probability and Statistics.: VIII. De Morgan and the Statistical Study of Literary Style.” Biometrika 45 (1/2): 282–82. https://doi.org/10.2307/2333072.
McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. Second. CRC Texts in Statistical Science. Boca Raton: Taylor and Francis, CRC Press.
McNeish, Daniel, Laura M. Stapleton, and Rebecca D. Silverman. 2017. “On the Unnecessary Ubiquity of Hierarchical Linear Modeling.” Psychological Methods 22 (1): 114–40. https://doi.org/10.1037/met0000078.
Mulvin, Dylan. 2021. Proxies: The Cultural Work of Standing in. Infrastructures Series. Cambridge, Massachusetts: The MIT Press.
Österman, Marcus. 2021. “Can We Trust Education for Fostering Trust? Quasi-experimental Evidence on the Effect of Education and Tracking on Social Trust.” Social Indicators Research 154 (1): 211–33. https://doi.org/10.1007/s11205-020-02529-y.
Senn, Stephen. 2003. “A Conversation with John Nelder.” Statistical Science 18 (1): 118–31. https://doi.org/10.1214/ss/1056397489.
Valentino, Nicholas A., Stuart N. Soroka, Shanto Iyengar, Toril Aalberg, Raymond Duch, Marta Fraile, Kyu S. Hahn, et al. 2019. “Economic and Cultural Drivers of Immigrant Support Worldwide.” British Journal of Political Science 49 (4): 1201–26. https://doi.org/10.1017/S000712341700031X.